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Abstract 

GPL Maxima is an open-source computer algebra system based on DOE-MACSYMA. GPL Maxima included 
two tensor manipulation packages from DOE-MACSYMA, but these were in various states of disrepair. One of 
the two packages, CTENSOR, implemented component-based tensor manipulation; the other, ITENSOR, treated 
tensor symbols as opaque, manipulating them based on their index properties. The present paper describes the 
state in which these packages were found, the steps that were needed to make the packages fully functional again, 
and the new functionality that was implemented to make them more versatile. A third package, ATENSOR, was 
also implemented; fully compatible with the identically named package in the commercial version of MACSYMA, 
ATENSOR implements abstract tensor algebras. 


1 Introduction 

GPL Maxima (GPL stands for the GNU Public License, the most widely used open source license construct) is the 
descendant of one of the world’s first comprehensive computer algebra systems (CAS), DOE-MACSYMA, developed 
by the United States Department of Energy in the 1960s and the 1970s. It is currently maintained by 18 volunteer 
developers, and can be obtained in source or object code form from http://maxima.sourceforge.net/ 

Like other computer algebra systems, Maxima has tensor manipulation capability. This capability was developed 
in the late 1970s. Documentation is scarce regarding these packages’ origins, but a select collection of e-mail messages 
by various authors survives, dating back to 1979-1982, when these packages were actively maintained at M.I.T. 

When this author first came across GPL Maxima, the tensor packages were effectively non-functional. The 
indicial tensor manipulation module could not display tensors, and could not carry out covariant differentiation; the 
component tensor manipulation package became nonresponsive the moment one tried to follow the documentation 
example and enter a metric. 

As it turned out, however, most of these problems were superficial in nature, and resulted from minor incompati¬ 
bilities between various LISP systems. In other words, they were easy to correct. Much of the work invested into the 
tensor packages since that time involves new and improved functionality. The purpose of this paper is to document 
the work that has been done. 


2 Tensors in curved spacetime 

Let us begin with a brief review of the tensor algebra topics that are implemented by the tensor algebra packages in 
Maxima. 

2.1 Tensor transformation laws 

The mathematics of curved spacetime can be said to be based on the transformation properties of various quantities 
under a change of coordinates (see [B] for an excellent introduction.) When a switch is made to a new coordinate 
system (indicated by primed coordinates below), the new and old components of a vector relate to each other as 
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follows: 


( X,1 \ 

x' 2 

\x' n ) 

Ordinary vectors that transform according to 0 are said to be contravariant. Loosely speaking, contravariant 
means “just like a vector”, and covariant means “just like a coordinate”. (By way of a simple example, consider 
a coordinate system that measures distances in meters. When you switch from meters to centimeters, your unit of 
measure decreases by a factor of 100. The quantities you measure, the length of a vector for instance, increase in 
contrast: a vector that was previously measured to be 5 units, i.e., 5 meters long, is now 500 units, i.e., 500 centimeters 
long. The coordinate units behave as covariant quantities, while the vector itself behaves as a contravariant quantity.) 

We can express the transformation laws for contravariant and covariant quantities in a more compact form. 
Denoting the transformation matrix in 0 by A, we have: 

x' = A ■ x (2) 

Now take another vectorial quantity, y, that is known to produce an invariant inner product with x: 

x-y = c (3) 

How would this product change if we were to switch from one coordinate system to another? Why, it’s supposed to 
be invariant, which means 

x' • y' = x • y = c (4) 

We know how to express x' as a function of x and the transformation matrix A. But the same transformation could 
not apply to y, because the result would not be the invariant quantity c: 

(A • x) ■ (A • y) ^ x • y (5) 

Then again, if you think about this a moment, the expression A • y doesn’t make much sense anyway. The vector 
x can be expressed in some coordinate system as a column vector. What do you multiply a column vector with in 
order to get a scalar? Why, a row vector of course. So y is a row vector. And how would you multiply a row vector 
on the left with the matrix A? Well, you could do so in principle, but the result is most certainly not another row 
vector. 

You can, on the other hand, multiply a row vector on the right with a matrix. And perhaps, instead of using A, 
you might consider using A , as in 

(A-x)-(y-A- 1 ) = xy (6) 

If this equation is satisfied for a quantity y, that quantity is said to be a covariant vector. A covariant vector maps 
a contravariant vector to the a scalar, namely the same scalar regardless of the coordinate system used. 

Many equations in physics and mathematics are matrix equations. For instance, we may have an equation in the 
form of 

v = T • u (7) 

When we change coordinate systems, we get: 

v' = A • v = A T u (8) 

How would we express v' as a function of u'? First, we multiply on the left by T _1 A _1 (remembering that matrix 
multiplication is not commutative but it is associative): 

T _1 ■ A^ 1 ■ v' = T" 1 • A -1 A T • u = T _1 T u = u (9) 

Next, we multiply on the left by A: 

A T 1 A -1 • v' = A u = u' (10) 
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Now that we managed to get u' on the right hand side, all that’s left is to eliminate all factors of v' on the left by 
multiplying both sides with A ■ T ■ A -1 : 

A • T ■ A- 1 • A • T- 1 ■ A’ 1 ■ v' = A ■ T ■ A- 1 • u' (11) 


Once again remembering that matrix multiplication associative, we can reduce this equation to: 

v' = (A • T • A- 1 ) • u' 

What we have here is a transformation rule for the matrix T: 

T' = A T A^ 1 


( 12 ) 


(13) 


We call a quantity like T a valence-) (j tensor. It has 1 contravariant index, and 1 covariant index. Generally, a 
quantity that transforms as 

X' = A fc • X • A~ l (14) 

is called a valence-[J] tensor. 

Since the matrix A -1 appears in these equations, it is a good question how this matrix can be computed. Inverting 
A the “hard way” would be an obvious method, but there is an easier way. Notice that A -1 expresses a vector 
transformation just as A does: 

x = A -1 • A • x = A -1 • x' (15) 

From this equation, we can read off the components of A -1 directly: 


ij = fed 


(16) 


2.2 Coordinate independence 

Many tensor operations exist that result in values that are independent of the coordinate system in which the tensors 
are expressed. We have already seen one example: the product of a covariant and a contravariant vector. If we 
specify the coordinate system, this product can be expressed as: 

n 

x-y = 'Y^x , y i (17) 

i=1 

where n is the dimensionality of the vectors, and we use upper indices to distinguish contravariant vectors from 
covariant vectors that have lower indices. 

The important thing here is that this product does not depend on our choice of a coordinate system. The result 
will be the same in any coordinate system. This observation is the basis of the so-called abstract index notation: 
expressing tensor equations using indexed quantities without explicitly specifying the coordinate system. 

Another powerful notational convention is the Einstein summation convenction. Basically, it means dropping the 
summation sign in equations in which an index appears in both a contravariant and a covariant position; in these 
cases, the summation is implied, i.e., 

n 

Y^x l y i =x l y i (18) 

This operation, in which a repeated tensor index (often called a dummy index) disappears as an implied summation 
is carried out, is called tensor contraction. 

2.3 The metric 

The distance between two points in a rectilinear coordinate system can be expressed using the theorem of Pythagoras: 

ds 2 = ^ (dx 1 ) 2 (19) 

i 
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When a generalized coordinate system is used in curved space, this formula obviously does not work. It does, however, 
apply so long as the quantities ds and dx 1 are infinitesimally small. How would the infinitesimal squared distance , 
ds, change under a change of coordinate systems? Here it is: 



l k i 

(20) 

By introducing the quantities 

dx n dx' 1 

9 kl 2^ Q x k g x i 

i 

(21) 

we can rewrite the formula as 

ds 2 = gudx k dx l 

(22) 

Notice that 

Qkidx = — t ~cjmidx dx = — 7 -ds 
ax K ax K 

(23) 


If we take the infinitesimal squared distance, ds 2 , to be unity, multiplying by gu accomplished a conversion from the 
infinitesimal quantity dx l to d/dx k . This is called a lowering of the index. What we basically have here is a unique 
way, with the help of g^i, to associate contravariant and covariant quantities. 

Similarly, the inverse of the metric tensor, g kl , can be used the raise a tensor index. The inner product of the 
metric tensor with itself gives the Kronecker-delta (i.e., the identity matrix): 

W m = 6 .r (24) 

The quantities gki themselves transform as a valence- [°] covariant tensor, which justifies the notation. 


2.4 Covariant differentiation 


In Euclidean space, when rectilinear coordinates are used, partial differentiation by a coordinate results in new 
tensorial quantities. For instance, partial differentiation of a scalar field by the coordinates defines the gradient of 
that field. Higher-order quantities can also be differentiated this way; we can, for instance, compute the gradient of 
a vector field and obtain a valence-tensor field. 

This is no longer true in curved space, where partial differentiation by a coordinate no longer necessarily produces 
a quantity that transforms as a tensor under a coordinate transformation. The question arises, then: is it possible to 
define a new type of a derivative operator, one that on the one hand has the same algebraic behavior (e.g., obeying 
the Leibnitz rule) as the ordinary derivative operator, but on the other hand, produces a quantity that behaves like 
a tensor under a change of coordinates? 

To answer this question, we first need to explore the concept of a parallel vector field. In flat space, a parallel 
vector field consists of identical vectors attached to each point in spacetime. In other words for a vector X 

vj 

dx ’ = Y.frX d * k =° (®) 


How does this expression transform under a change of coordinates? We can answer this question by utilizing the 
identity X' j = Y^ h | &rX h , to get: 


h k 


(26) 


With the help of the metric, it is possible to eliminate the second derivatives from this equation. First, we need to 
differentiate eu with respect to x h and then cyclically permute the indices to obtain the following identities: 

dgkl _ V ( dx ,j , dx ,j d 2 x' d \ 

dx h 2-*jj \ dx h dx k dx L ' dx k dx h dx L J 

dgih ( d 2 x ,j dx'i , dx'i d 2 x /J \ 

dx k 2-^dj y dx k dx l dx h ' dx 1 dx k dx h J 

dffhk — ST' ( d 2 x ,j dx ,j I dx ,j d 2 x' J \ 

dx 1 Z—sj ydx l dx h dx k ' dx h dx l dx k J ' ' 
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Adding the first and second of these three equations, subtracting the third, and dividing the result by two, we 
obtain: 

_ 1 / dg k i dgih dg hk \ _ <9 dx ,J 

2 \ dx h dx k dx l ) dx h dx k dx l 


hkl 


(28) 


Now we multiply by dx'i/dx l and sum over j to get 


dx' j 


E g x i dX J E ghidX + 


(29) 


h k 


Both sides of this equation should be zero for a parallel vector field. The right hand side can be made more meaningful 
if we multiply it by g^ 1 (we’re also dropping the summation signs now, relying on the summation convention instead): 

dX j + g jl T hk iX h dx k = 0 (30) 

In other words, for a parallel vector field, the quantities Fhk J = g^Fhki express the difference by which the ordinary 
differential operator will fail. This suggests a definition for a differential operator in curved space: 

BXi 

V fc X? = ~dx k+Thk3Xh (31) 

For a covariant vector field, the same derivation yields 

8X 

XkX j = -g^- r jk h X h (32) 

And, for an arbitrary tensor of valence- A J, we get 




i 1 ■■■l, 

dx k 


E r - 
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E r. mrpji...j r 

L l pk J -i 1 ...ip_ 1 mi l 3 + 1 ...i a 


(33) 


The quantities T^-fc and T'ij k are called the Christoffel-symbols of the first and second kind, respectively. 

The operator we just defined, V&, can justifiably called a differential operator for several reasons. First, as can 
be verified by direct computation, it obeys the Leibnitz rule: V k (X • Y) = (VfcX) • Y + X • (V*,Y). Second, it can 
be shown that once the metric, gij, is given, the differential operator is uniquely defined, in that it is the only such 
operator that preserves the inner product of two vectors parallel transported along any curve, i.e., XkgijV l w J = 0. 

Note that a shorthand is often used to express the covariant derivative; for instance, instead of V k Xij, we often 
write Xij - k . 


2.5 The Lie-derivative 


The covariant derivative is tensorial in nature, but unfortunately, it depends on the metric. The ordinary derivative 
is not dependent on the metric, but it’s not tensorial: it depends on the choice of the coordinate system. So the 
question arises: is it possible to construct operators that are not dependent on the metric, but produce tensorial 
results that are not coordinate system dependent? 

There are, in fact, many such operators possible (see [7]), most notably among them the Lie-derivative: 


C v Ht;± = V h d h Ht k - HfcfdhV* - ... - Hl;±d h V k 


Hj;±d l V h 


H\;; k d n V h 


(34) 


The geometric meaning of the Lie-derivative is directional differentiation: it tells you how a tensor field changes in 
the direction of the vector field V. 


2.6 Commutators 

The ordinary partial differential operator is known to commute: dj,d 7 X — djd, X = 0. Is this also true for the covariant 
derivative? 
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For a scalar field /, the commutator is 




d 2 f r » df d 2 f 

dx k dx l lk dx h dx k dx l 


V h ^L - (r h -r M df 

kl dx h ' U lk ' dx h 


(35) 


The way we defined them, the Christoffcl-symbols are symmetric in their first two indices, and therefore T k i h — Ti k h = 
0. We could conceivably define another differential operator whose connection coefficients (the generalization of the 
concept embodied in the Christoffel-symbols) are not symmetrical, and thus c k i h — ci k h does not vanish; this quantity 
is then called torsion. 

How about the commutator of the covariant derivative on a vector field? Here it is: 


(V fc V,-VzV fc )* J ‘ = 

(VW J ) + iW (V*x m ) - T lk m (V m P) - A (VfcJW) - T mt j ( V k X m ) + (V m X J ) = 

(aSw - ») + a? ( rv 'W - 5? + 


k J 
dr h f 


'x h + (r kl m -r lk m )(v m x j ) = 


hr x h + T hl jdx 


dri 


hk 


dx k 

mk^ r hi 

dT hl i dT hk 3 
dx k 


x n - r 


, dX h 


U I ^ 1 ■*" hk r\ 1 1 TTLK, r\ 1 

OX K OX 1 OX 1 OX L 

P jp m \rh p 7p m \rh , /p m p m\ (\ 7 vi\ 

J- rrt.kf 1 hi A -A +(lfcZ — 1/fe j (V m A J J = 




9a: fe 


9a; z 


+ iwr w m - r mr T hfc m + (r fei m - r ifc m ) (v ro ?w) 


(36) 


The second part of this result we already recognize: it’s the torsion again. The first part contains the Riemann 
curvature tensor, defined as 


Rhik 1 = ~Rhkl 3 = -+ T mk J T h i m - T ml >r hk m (37) 


Several other quantities are derived from the curvature tensor that are widely used in physics and mathematics. 
These include the Ricci-tensor: 



Rij — Rji — Rijk 

(38) 

the scalar curvature: 


R = R\ 

(39) 

and the Weyl-tensor CpTil: 

W^ijkl — Rijkl ^ 

2 2 

l)(n — n _ 2 9j[iRi]k) 

(40) 


(Indices in square brackets indicate complete antisymmetrization; for instance, Huy = (H^ — Aji)/2\.) 

In gravitational theory, the Ricci-tensor expresses the contribution of matter and energy to the curvature tensor, 
while the Weyl tensor (often called the conformal tensor, and labelled C a b c d ) represents the curvature of empty 
spacetime, e.g., curvature far from any matter and energy sources, caused by free gravitational waves. 

Yet another important tensor in gravitational theory is the Einstein tensor: 

Gij = Rij - i Rg i:j (41) 


2.7 Torsion and nonmetricity 

As described in the previous section, if the connection coefficients are not symmetric in their indices, the metric is 
said to have torsion. 
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Torsion can be characterized by the torsion tensor T-ij k . When torsion is present, it modifies the connection 
coefficients such that: 

t'l-jk = T ijk — (Tjj (Jkm T Tki Qjm T Tfcj 9'i.rn ) (42) 

The quantites used to modify the Christoffel symbols in this expression are often called contortion coefficients: 


k 'i'jk 2 (Ay 9km + T k i m gjm + T kj m g im ) 


(43) 


Using the derivative operators based on the Christoffel-symbols (with or without torsion) we find that the covariant 
derivative of the metric tensor is identically zero: 


V fc5lj = 0 (44) 

It is possible to modify the derivative operator, while retaining its algebraic properties, so that the covariant derivative 
of the metric tensor is no longer zero. This property of a derivative operator is called nonmetricity. Nonmetricity is 
characterized by the vector field /z/.: 

gij^k = gkgij ( 45 ) 

The contribution of the nonmetricity field to the connection coefficients can be expressed as 

C-ijk — h ijk T — (gikjlj 4“ gjkj^i gijJ^k) (45) 

The quantities that appear in this equation are sometimes called nonmetricity coefficients: 

k'ijk “ ( 9ikRj gjkt^i T 9ij P k) (47) 


2.8 Rigid frames 

In a smooth space, the immediate vicinity of every point can be approximated by a flat space (i.e., a space without 
curvature.) This gives rise to the idea of attaching a rigid (typically, though not necessarily, rectilinear) set of basis 
vectors (a rigid frame , also called (in four dimensions) a tetrad, vierbein , or, in arbitrary dimensions, a vielbein ) to 
each point. While this can be done at each point, the resulting collection of basis vectors will not form a coordinate 
system. Intuitively, this means that supposedly parallel lines in a coordinate grid will, sooner or later, meet and 
cross each other. Such a basis is said to be nonholonomic. 

Though not a coordinate base, such a base of moving frames can nevertheless be used to measure vectors and 
other quantities in that space. And, at least in some cases, the use of a nonholonomic frame can result in simpler 
equations. 

A nonholonomic base is characterized by a set of n basis vectors e( Q ), where n is the dimensionality of the space, 
and parenthesized indices refer to the basis vectors themselves. Each of these basis vectors can have (contravariant) 
components relative to some coordinate base: e ( a y The inner product of basis vectors defines the frame metric: 
e (a) e (b)i = Vab- Typically, the frame metric would either be the unit matrix (for a Euclidean frame) or a unit matrix 
with Minkowski signature (for a Lorentz frame.) 

Raising and lowering of the frame basis indices can be accomplished using the frame metric or its inverse; e.g., 
e (a) _ rj ab e ^ b y Since g a bg bc = it’s also true that = 5^. 

To compute the covariant derivative of a quantity expressed using these basis vectors, we need the Ricci rotation 
coefficients (0): 

'tabc = e (a)i;k e (b) e (c) (48) 

In a frame base, these coefficients take the place of the Christoffel symbols. Computing them in this form is 
not practical, however, since we’d need a way to compute the covariant derivative of the basis vectors before the 
coefficients are calculated. What we do instead is take a linear combination of these coefficients: 

^abc = ^)abc 'Yacb = (^(a)i;fc ^(a)fc;z) ^(b)^(c) (49) 

The quantities \ a bc or A 0 b c are often called the frame bracket. 
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The last expression can be further simplified, with ordinary partial differentiation replacing the covariant deriva¬ 
tive, as the Christoffel symbols cancel out due to their index symmetries, leaving only the torsion: 

(^( a)i;k ^(a)fc;z) ^(b)^(c) (^(a)z,fc ^(a)k,i ^ik &(a)m) (50) 

Using these quantities, the Ricci rotation coefficients can be expressed as 

'Yabc = 2 abc T ^bca ^ cab ) (51) 

The Ricci rotation coefficients are antisymmetric in their first two indices: 'yabc = —7 bac- This means that in contrast 
with the n 2 {n + l)/2 independent components of the Christoffel-symbols, the Ricci rotation coefficients have only 
n 2 (n — l)/2 independent components: 24 instead of 40 in the case of 4 dimensions. This explains why sometimes, 
the use of a frame base presents significant advantages. 

Expressed using the Ricci rotation coefficients, the Riemann tensor takes the following form (cf. TOJ, eq. 3.4.21): 

Rdabc — ^ (a)^ i^fbcd C^Vj7acd "t“ TJ ^ (7 a fobbed 7fr/c7aed T ^/afb^/ecd 7fr/a7ecd) (52) 


2.9 Algebraic classification of metrics 

When studying curved manifolds, one frequently encountered problem is whether or not two metrics represent the 
same manifold. Usually it is not at all obvious from the algebraic properties alone if two metrics are equivalent; 
consider, for instance, the metric associated with polar vs. rectangular coordinates ( dx 2 + dy 2 vs. dr 2 + r 2 d(j) 2 ) or, 
in 3 dimensions, Cartesian vs. spherical coordinates ( dx 2 + dy 2 + dz 2 vs. dr 2 + r 2 d9 2 + r 2 sin 2 9d(j) 2 ). A variety of 
techniques exists to help establish if two metrics are guaranteed to be different or may be the same; most well known 
among these, perhaps, is the Petrov-classification, applicable to 4-dimensional metrics with Minkowski signature. 

The Petrov-classification is based on the fact that the Weyl-tensor of any metric gives rise to several invariant 
quantities. Denoted by the symbols 9/ 0 , ... T 4 , these invariants provide the coefficients of a quartic equation for 

the principal null directions of the metric (see 0 for details.) 

To compute the Petrov classification, one requires two elements: the Weyl tensor in a coordinate base, and a 
Newman-Penrose null tetrad representing the metric. For instance, if we have an orthonormal tetrad base e( a ) with 
(+, —, —, —) signature and coordinates (f, x, y, z), the corresponding null tetrad can be computed as: 

k = ^ (e 1 + 62 ) 
l = ^{e 1 - e 2 ) 
to = ^ (e3 + *64) 

to = -^ (e 3 — TO 4 ) (53) 


With the Weyl-invariants and the null tetrad at hand, classification of the metric can proceed. Two principal 
quantities used in this computation are 

I = - 4T4T3 + 3^2 (54) 


and 


J = 


T 0 

Ti 

T 2 


Ti 

T 2 

^3 


T 2 

^3 

T 4 


(55) 


The actual algorithm is quite messy: the Petrov type is determined by which of the T invariants are non-zero; in 
some cases, this pattern directly determines the Petrov type, while in other cases, /, J, and additional coefficients 
need to be computed first. A complete description of a modern algorithm can be found in (|HJ)• 

The algebraic types obtained by the Petrov-classification can be represented in a diagrammatic form illustrating 
increased specialization: 


I 




These algebraic types also possess geometric significance. They determine the number of principal null direc¬ 
tions that a metric possesses. A general four-dimensional spacetime has four principal null directions, but in the 
algebraically special cases, some of these null directions coincide: 


Type 


Null directions 

I 

(1,1,1,1) 

All null directions are independent 

D 

(2,2) 

Two pairs of coinciding null directions 

II 

(2,1,1) 

Two null directions coincide, two are independent 

III 

(3,1) 

One triplet of coinciding null directions, one independent 

N 

(4) 

All null directions coincide 


For a more detailed explanation, see (7 and Sj. 


2.10 Exterior forms 

A particularly interesting variety of tensors are fully antisymmetric covariant tensors. I.e., tensors of the type 

A i...k...l...n A 

These tensors are notable because they form the basis of Cartan’s exterior calculus of differential forms. Simple 
operations exist on these tensors that can be expressed in an index-free notation. They include the exterior derivative, 
contraction with a vector, and the exterior, or wedge product. 

The wedge product is the fully antisymmetrized product of two covariant, fully antisymmetric tensors. Unfor¬ 
tunately, not all authors agree on the precise definition of the wedge product. As Penrose and Rindler comment 
in “We should remark that, in the literature, a slightly different convention is frequently employed, in that the 
quantities a ai ... ap := p!A Ql ... ap rather than our A ai ... Q are used to denote the components of a p-form.” 

The advantage of the notation employed by Penrose and Rindler, and also used in El and J3, is that it agrees 
with the concept of total antisymmetrization. The other notation, used for instance in [3| or J2] , is more geometrically 
inspired, as it produces exterior forms that correspond with the notion of the volume element. 

Generally speaking, in works about tensor algebra, the definition employed by Penrose and Rindler is used more 
frequently; in works utilizing the index-free notation of exterior forms, the “geometric” notation is preferred. 

How does this affect the wedge product? If the “tensorial” notation is employed, the wedge product of a p-form 
and a (/-form is defined as: 


A ii...i p A Bj 


1 


.. .kpll.. .lq 


In the “geometric” case, however, the definition is altered: 


kl...kn 


Bu 


A n...u. A : A *...k P B h...i q 


(56) 


(57) 


The exterior derivative is formally defined as the wedge product of the partial derivative operator and a tensor. 
Consequently, this definition is also dependent on how we define the wedge product. In the “tensorial” case, the 
exterior derivative is defined as 


dA — ^[i A j 1 ...j p ] 


1 

(j>+ 1 )! 


rkl i.. .Ip 


V k Ajj... i p 


(58) 


In the “geometric” interpretation, the definition of the exterior derivative changes to 


dA — (p+ 1 )V[iAj 1 ...j p ] 




(59) 


Fully antisymmetric covariant tensors, like any other tensors, can be contracted with a vector: e.g., v What 
makes the case of fully antisymmetric covariant tensors unique is that the index by which the contraction takes place 
does not matter: if a different index is used, only the sign of the result will be affected, and the sign change will be 
consistent with the permutation of indices. 
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2.11 Abstract tensor algebras 

Tensors, and matrices representing tensors, offer rich algebraic structures. As a simple example, 2 x 2-nratrices in 
the form “) can be used to represent the algebra of complex numbers. Other examples include the Pauli and 
Dirac matrices, or the algebra of quaternions. 

One of the simplest algebra types is the symmetric algebra , defined by the commutation rule u ■ v — v ■ u = 0. 

Grassmann algebras are anticommutative: u ■ v + v ■ u = 0. The direct sum of all exterior forms in dimension n 
forms a Grassmann algebra. 

A Clifford algebra (see [J) is defined by the scalar unit, one or more basis vectors, and a symmetric scalar 
anticommutator: u ■ v + v ■ u = 2f s (u,v). Perhaps the simplest example is R(0,1), a Clifford algebra with 1 basis 
vector (iq), and an anticommutator function that is defined as f s (v i,iq) = — 1. This is none other than the algebra 
of complex numbers. 

Similarly, the Clifford algebra R(0, 2) with a anticommutator defined as f s (vi,vi) = f s (v 2 ,^ 2 ) = — 1, and 0 for 
all other arguments, defines the algebra of quaternions; V\, V 2 , and their product, V\ ■ V 2 will correspond with the 
three quaternionic imaginary units. 

Other, notable Clifford algebras include the algebra of Pauli spin matrices (R(3,0)) and Dirac spin matrices 
(R(3 ’ 1)} - 

Closely related are the symplectic algebras : these are defined by an antisymmetric scalar commutator function, 
such that u ■ v — v ■ u = 2 f a (u, v). 

Lie enveloping algebras are characterized by the Jacobi-identity: [it, [iquj] + [it, [k;, it]] + [w, [it, i>]] = 0, where [, ] 
denotes the commutator, an antisymmetric vector-valued function: [u,v] = u.v — v.u = 2v a (it,n). 


3 Tensors in Maxima 

A computer algebra system can represent tensor quantities in many ways. 

First, tensors can be represented in terms of their components relative to some coordinate system. Second, 
tensors can be represented by their index properties. Third, tensors can be represented by symbols, with special 
simplification rules used to resolve expressions containing these symbols. 

All these representations are implemented by the various Maxima tensor packages: ATENSOR, CTENSOR, and 
ITENSOR. 

3.1 Indicial tensor manipulation 

The ITENSOR package treats tensors as “opaque” objects. Tensor values do not matter; tensor expressions are 
evaluated using formal algebraic rules, in particular algebraic rules concerning abstract tensor indices. In order to 
understand how this package works, it is nececssary first to take a brief look at how Maxima represents algebraic 
objects. 

3.1.1 The Maxima architecture 

It is not the purpose of this paper to provide a comprehensive overview of the Maxima system. (The author must 
admit that many areas of this large, complex system remain to him a complete mystery.) What is relevant with 
respect to the ITENSOR package is how Maxima represents expressions, function calls and index expressions in 
particular, using the LISP language. 

The core concept of the LISP language is, unsurprisingly, the idea of a list. Or, to be more precise, the idea of 
an ordered pair, the first element of which is the head (or car), the second element the tail (or cdr) of the list. List 
elements are themselves either lists or atoms: e.g., a number, a symbol, the empty list (nil). 

For instance, the list (12 3) would be represented as the ordered pair of the elements 1, and (2 3). The first 
of these is an atom, the second itself a list, i.e., another ordered pair. The first element of this second pair is the 
atom 2, the second the list (3). This last list, too, is represented by an ordered pair: the first element is the atom 
3, the second element the empty list nil. In diagram form: 
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(1 2 3) 


1 




2 


(2 3) 

\ 

(3) 

/ \ 

3 nil 


LISP provides operator functions that are used to extract elements or sublists from a list. The most common of 
these are car (to obtain the head of a list) and cdr (to obtain the tail). In other words, car gives you the first of 
the ordered pair of elements, cdr the second. 

All LISP implementations provide many convenience functions that are combinations of car and cdr. For instance, 
cadr is the head of the tail of a list: e.g., (cadr ’ (1 2 3)) gives the atom 2. In other words, (cadr x) is shorthand 
for (car (cdr x)). 

LISP is a particularly suitable choice of a programming language for representing mathematical expressions. LISP 
has a strict “prefix” syntax: operators, be they unary, binary operators or function calls, always precede their list of 
arguments. As an example, here is a possible LISP representation of the mathematical expression x + ysinz: 

(+ x (* y (sin z))) 

Maxima uses similarly constructed LISP lists to represent expressions. However, instead of LISP operators, which 
may not always perform in the desired manner, Maxima implements its own set of operator functions. So the same 
mathematical expression looks slightly different as a Maxima object: 

((mplus) $x ((mtimes) $y (($sin) $z))) 

The dollar sign in these expressions simply indicates a Maxima variable, as distinguished from LISP-level keywords 
or variables. 

As a computer algebra system, Maxima has a powerful simplification capability. When an expression is entered 
by the user, it is automatically processed by the simplifier. To avoid unnecessary processing, the simplifier “marks” 
expressions in various ways. After being processed by the simplifier, the previous expression will appear thus: 

((mplus simp) $x ((mtimes simp) $y (($sin simp) $z))) 

Maxima also has lists. These are to be distinguished from LISP lists. In the syntax of the Maxima language, a 
Maxima list is a comma-separated list of objects enclosed in square brackets. Internally, the Maxima list [a,b,c] is 
represented in LISP as: 

((mlist) $a $b $c) 

3.1.2 Indexed objects 

Tensors are physical/geometric quantities that are characterized by their transformation rules under coordinate 
transformations. Though it is often convenient to think of tensors as matrices, a matrix is nothing more but the 
representation of a tensor in a particular coordinate system; when the coordinate system changes, so do the elements 
of the matrix, even though the tensorial quantity they represent remains the same. 

Nevertheless, a number of tensorial operations exist that can be carried out in an arbitrary coordinate system, 
yet the result would be a proper representation of the appropriate tensorial quantity in that coordinate system. This 
means that a coordinate transformation can be applied either before or after the operation is performed, and the 
result would be the same. Symbolically, if / is a tensor operator acting on a set of tensors T), T 2 ... T„, and ' denotes 
the coordinate transformation, the following is sometimes true: 

\j{T 1 ,T 2 ...T n )]' = f{r 1 ,r 2 ...r n ) (eo) 

In particular, / can be tensor addition, the inner or outer product of tensors, and certain types of tensor differenti¬ 
ation. 


11 



This gives rise to the abstract index notation which one frequently encounters in the literature: we do not care 
what coordinate system is used, but we employ indices to unambiguously specify the operation. 

Abstract indices can be covariant and contravariant. 

One reason why the abstract index notation proves to be a very powerful algebraic tool is a “feature” of the 
inner product operation: index contraction. For instance, it is known that the metric tensor, g, contracts with other 
tensors, and can be used to raise or lower tensor indices. So for instance, the following is treated as an identity: 

g ab T bc = T c a 

Notable is the fact that this operation can be carried out without knowing anything about what’s “inside” the 
tensor T: the rules are strictly formal in nature, for manipulating the “indexed object” T, and the only information 
needed to carry out the operation is the pattern of indices that the object possesses. 

3.1.3 Representing tensors in ITENSOR 

It is this formalism that is captured by the ITENSOR package. At the risk of sounding like a marketing brochure, it 
must be mentioned that it is this capability that makes ITENSOR particularly unique and powerful: it can be used 
for many tensorial problems very efficiently, as no computational resources are wasted to manage the “contents” of 
a tensor when, in fact, algebraic manipulation of the indexed symbols can suffice. 

The ITENSOR package represents tensors as formal Maxima function calls. For instance, the tensor T£ would 
be represented as: 

T([a] , [c]) 

Many of the features of ITENSOR can be discerned from this simple example. The tensor T is formally represented 
as a function call with two parameters. Both parameters are Maxima lists: the first contains the covariant, the second, 
the contravariant indices of the tensor. 

The ITENSOR package includes a function whose purpose is to visually present such indexed objects in the 
familiar notation. This is illustrated by the following excerpt from a Maxima session: 

(7,i2) ishow(g([a,b] , []))$ 

(“/.t 2) g 

a b 

An experienced user would frequently use the ishow function inside tensorial 
interim results in a visually pleasing form. 

As remarked earlier, a Maxima function call is represented internally as a LISP 
is represented internally by ITENSOR as follows: 

((mlist) $g ((mlist) $a $b) ((mlist))) 

Note that although the list of contravariant indices is empty, it is nevertheless present. 

Tensors in ITENSOR can also have derivative indices. These indices represent ordinary partial derivatives in the 
covariant position. For instance, the tensor g a b, c would be expressed using ITENSOR as g( [a,b] , [] , c). (Note that 
the similar, widely used notation to represent covariant differentiation, e.g., g a b;c, is not supported by the current 
version of ITENSOR.) 

The LISP list corresponding with the tensor g a b,c is 
(($g) ((mlist) $a $b) ((mlist)) $c) 

This formalism makes it particularly easy to extract the list of covariant, contravariant, and derivative indices of 
a tensor. Once it has been ascertained that, say, $x is a tensor object, the LISP expressions (cdadr $x), (cdaddr 
$x), and (cdddr $x) can be used to obtain the list of covariant, contravariant, and derivative indices, respectively. 
Tensor contraction is demonstrated by the following excerpt: 


expressions in order to present 
list. Accordingly, the tensor g a b 
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done 


(7.13) imetric (g) ; 

(7.03) 

(7.14) isfiow(g([b,c] , [])*T([], [a,b]))$ 

a b 

(7.t4) T g 

b c 

(7.15) ishow (contract (7.) ) $ 

a 

(7.t5) T 

c 

Incidentally, this example already demonstrates one of the shortcomings of the notation employed in ITENSOR. 
When the index just lowered is raised again, we should get back the original tensor: i.e., g dc gb c T ab = g dc T a c = T ad . 
However, when we carry out the computation in Maxima, we get a different result: 

(7.16) ishow(contract (g( [] , [d, c] ) *7.) ) $ 

d a 

(7.t6) T 

The reason for this discrepancy is due to the fact that ITENSOR maintains contravariant and covariant indices 
in two separate lists. When a contravariant index is lowered, information regarding its position relative to other 
contravariant indices is lost; when the index is subsequently raised, its position among contravariant indices will be 
predetermined by the contraction algorithm, and will not in any way be related to its location prior to its lowering. 

3.1.4 A new index notation 

This is one of the shortcomings in ITENSOR that was addressed by the present work. The goal was to develop a 
notation that is a) fully compatible with the current ITENSOR notation and can be used interchangeably with it; 
b) is relatively easy to implement, requiring only minor modifications to ITENSOR’s complex algorithms; and c) 
provides a tensor notation that can be used to keep track of index ordering in arbitrarily complex tensor expressions. 

Sometimes, it is best not to reinvent the wheel. For this reason, the present author examined how this problem has 
been addressed in other tensor algebra packages. For instance, the built-in tensor algebra package in Maple 9.5 uses 
the relatively cumbersome method of explicitly spelling out the position of indices using the “index characteristic” 
component of a tensor object. Tensor objects are thus defined using expressions like this one: 

> g:=create([-1,-1],some_array); 

Though unambiguous, this notation does not make it easy to enter and manipulate complex equations involving 
many indexed objects. 

The downloadable grTensorll package for Maple, designed specifically for use in general relativity, employs a 
similar notation. 

The Derive package, currently marketed by Texas Instruments, uses a somewhat more efficient representation. 
In Derive, the underscore character in the tensor’s name is used to separate groups of covariant and contravariant 
indices. For instance, the covariant metric tensor, g.y, in Derive is denoted as g_ij, while the (1,3) Riemann tensor 
R l jki would be denoted as R__i_jkl. 

Both packages suffer from one serious shortcoming, however. They do not conceptually separate the components 
of a tensor from the algebraic symbol of the tensor. As has been remarked above, many tensorial equations can be 
computed utilizing only the algebraic properties of the tensors involved; carrying the components, in these cases, 
would mean an unnecessary computational burden. 

In any case, while each of these notational conventions has merits, they are not directly applicable to ITENSOR, 
since they substantially differ from the indexed object notation employed by this package. It is obvious that whatever 
convention is used, the basic formalism of a tensor as a formal function call with indices supplied in the form of list 
arguments should be preserved. 

As such, two possible extensions to the tensor syntax present themselves. One possibility, inspired by Derive, is 
to allow multiple groups of covariant and contravariant indices. For instance, the tensor T a b c d could be represented 
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as the Maxima object T ([a] , [b] , [c] , [d]) where each list in an odd position contains covariant indices, while lists 
in even positions contain contravariant indices. Index raising and lowering would move indices only to neighboring 
lists, which would ensure that an index that is raised and then subsequently lowered retains its original position 
relative to other indices. 

The main shortcoming of this notation is that it is still relatively cumbersome to use. It also makes it somewhat 
difficult to extract derivative indices; the algorithm would have to enumerate all parameters to the formal function 
call that represents a tensor, find the first parameter that is not a list, and extract this and the remaining parameters. 
It would also be quite difficult to adapt many of the utility functions of ITENSOR to successfully deal with this 
notation. 

For this reason, eventually another notation was chosen. To understand the rationale for this choice, first notice 
that in the abstract index formalism indices are truly abstract: they merely serve as labels, nothing more. In 
particular, no valid tensor equation contains mathematical expressions in the place of an index: T l+ \ is not a valid 
symbol. 

What this means is that we are free to use some mathematical operator to “mark” indices. Specifically, nothing 
prevents us from using the negative sign, for instance, to mark indices that are in the “wrong” list: e.g., contravariant 
indices that appear in the list of covariant indices, or vice versa. 

In other words, in this notation the tensor T a b c d could be represented as the object T( [a,-b,c,-d] , []). This 
would be effectively equivalent to the “old” notation of T([a,c] , [b,d]) with one crucial difference: the relative 
ordering between covariant and contravariant indices is preserved. 

How difficult is this notation to implement, while permitting it to coexist with the “old” notation? Fairly easy, as 
it turns out. For most ITENSOR functions, the actual index ordering doesn’t matter, so in these cases, all we need is 
new functions to extract the list of covariant and contravariant indices. These tasks can be accomplished easily using 
two helper functions that extract the “positive” and “negative” indices from a list. In order to understand these 
helper functions, it is necessary to know how Maxima represents the negative sign: what is used is the expression 
— 1 • x, as in ((mtimes) -1 $x), to represent —x. Thus we have the following LISP code: 


(defun plusi(l) 

(cond 

((null 1) 1) 

((atom (car 1)) (cons (car 1) (plusi (cdr 1)))) 

((and (eq (caaar 1) ’mtimes) (eq (cadar 1) -1)) (plusi (cdr 1))) 
(t (cons (car 1) (plusi (cdr 1)))) 

) 

) 


(defun minusi(l) 

(cond 

((null 1) 1) 

((atom (car 1)) (minusi (cdr 1))) 

( 

(and (eq (caaar 1) ’mtimes) (eq (cadar 1) -1)) 
(cons (caddar 1) (minusi (cdr 1))) 

) 

(t (minusi (cdr 1))) 

) 

) 


With these two functions, obtaining the list of covariant and contravariant indices is extremely simple: 

(defun covi (rp) (plusi (cdadr rp))) 

(defun conti (rp) (append (minusi (cdadr rp)) (cdaddr rp))) 

Armed with these two functions, all that was left to be done is to substitute (covi x) and (conti x) in place of 
(cdadr x) and (cdaddr x), respectively, in all the appropriate places. (A slight complication arises from the fact 
that in some cases, the “pure” list of indices was obtained, while in other cases, a Maxima-style list which begins 
with the object (mlist) was needed, but this is an implementation detail.) 
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The place where index order does matter is inside the contract function. More specifically, inside the (internal) 
helper function contract 1. The following diagram depicts the various contraction functions that implement index 
contraction for different object types: 

$C0NTRACT 

I 

I —C0NTRACT5 

I 

+—C0NTRACT4 

I 

+-C0NTRACT3 

I I 

I +-C0NTRACT1 

I I 

I +-C0NTRACT2 

I 

+-C0NTRACT5 

As this diagram shows, the actual contraction of two tensors is carried out by contract 1 . Here is the relevant 
fragment of this function that checks if there are covariant indices (c) in an object that match the contravariant 
indices (b) of the contracting object, or contravariant indices (d) in an object that match the covariant indices (a) 
in the contracting object: 

(cond 

( 

(and b c (setq f (contracts b c))) 

(setq b (car f) c (cdr f)) 

) 

( 

(and a d (setq f (contract2 ad))) 

(setq a (car f) d (cdr f)) 

) 

(The function contract2 takes two lists of indices as its parameters and returns a pair of lists: the car contains 
the “free” indices in the first argument, while the cdr contains the “free” indices in the second argument.) 

This area of the code needed to be expanded, in order to accommodate the new style indices. Two separate cases 
had to be addressed: when the covariant list of the contracting object contained any “negative” indices that matched 
covariant indices in the other object, or if the covariant list of the contracting object contained any “positive” indices 
that appeared as “negative” indices in the covariant list of the other object. These two tasks are accomplished by 
the following code fragment: 

(cond 

( 

(and a (minusi c) (setq f (contract2 a (minusi c)))) 

(do 

( 

(i c (cdr i)) 

(j (car f)) 

(k) 

) 

((null i) (setq a (removenotin j a) c (reverse k))) 

(cond 

( 

(or (atom (car i)) (member (caddar i) (cdr f))) 


Called for RPOBJ’s; contract a tensor with itself 

Called for product expressions 

Contract an object with an element of a list 

Contract two tensors 

Remove like members from two lists 
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(setq k (cons (car i) k)) 


) 

( 

(not (null j)) 

(setq k (cons (car j) k) j (cdr j)) 

) 

) 

) 

) 

( 

(and (minusi a) c (setq f (contract2 (minusi a) c))) 

(do 

( 

(i c (cdr i)) 

(j (car f)) 

(k) 

) 

((null i) 

(setq 

c (reverse k) 
a (append 
(plusi a) 

(mapcar #’(lambda (x) (list ’(mtimes simp) -1 x)) j) 

) 

) 

) 

(cond 

((member (car i) (cdr f)) (setq k (cons (car i) k))) 

( 

(not (null j)) 

(setq k (cons (list ’(mtimes simp) -1 (car j)) k) j (cdr j)) 

) 

) 

) 

) 

(As its name suggests, the removenotin (i 1) helper function removes objects not in the list i from the list 1 
and returns the truncated list as a result.) 

Surprisingly, one area of the ITENSOR package did not need to be modified: simplification functions that use 
the symmetry properties of tensorial objects to simplify tensor expressions. 

In ITENSOR, it is possible to declare a tensor to be symmetric in some of its indices using the decsym function. 
For instance, the metric tensor is typically symmetric in both its covariant and in its contravariant form. Assuming 
that g is the metric tensor, these symmetry properties would be expressed using the following commands: 

decsym(g,2,0,[sym(all)],[]); 
decsym(g,0,2,[],[sym(all)]); 

This notation remains fully applicable when “new” style tensor objects are used. In fact, the new notation makes 
things simpler for the Maxima user. The single declaration 

decsym(g,2,0,[sym(all)],[]); 

is sufficient to express that g a b = gba, g ab = g ba , and g a b = g b a , so long as all indices are placed in the first, 
“covariant” list, with a minus sign used to denote contravariant indices. 
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3.1.5 Riemannian geometry 

The ITENSOR package contains built-in definitions for the Christoffel symbols of both kinds, and the (1,3) Riemann- 
tensor, expressed in terms of the metric tensor. 

With the help of the Christoffel-symbols, ITENSOR could express the covariant derivative of any tensor field, 
using the covdiff function. This function, unfortunately, was broken in ITENSOR. The work to make ITENSOR 
fully functional again began here, by fixing this subroutine. 

The first new feature that was implemented in this area was the ability of ITENSOR to handle differentiated 
forms of the Christoffel-symbols. Previously, an expression in the form of, say, ichr2( [a,b] , [c] ,d) resulted in an 
error. Now, however, these forms are correctly evaluated in terms of the metric tensor. 

A far more ambitious set of changes was implemented on the basis of the reference manual of the commercial 
version of MACSYMA (0]). The ITENSOR package in this product possessed the ability to perform computations 
not only in a coordinate base but also in a frame base; furthermore, it had the ability to deal with torsion and 
nonmetricity. 

A frame base is characterized by two entities: the basis vectors e( Q ^, which can be collectively expressed in 
tensorial form, and the frame metric r/ a b- The latter would normally be either the identity matrix (for a Euclidean 
orthonormal frame) or the Lorentz-matrix (for a frame with Minkowski signature.) 

Note that ITENSOR does not have the capability to distinguish frame and coordinate indices. This can lead to 
ambiguities in ITENSOR expressions. This represents few problems so long as 1) the same symbol is not reused 
as both a frame and a coordinate index, to avoid contraction between incompatible indices, and 2 ) differentiation 
is carried out on expressions containing either only coordinate or only frame indices, with the iframejflag flag set 
appropriately. See the discussion at the end of section 2.3 of 0 for details. 

When a frame base, torsion, or nonmetricity is used, in most expressions involving curvature, the Christoffel- 
symbols are replaced by connection coefficients: 

Cabc — Tabc ^afec ^a£>c ( 61 ) 

when a coordinate base used, or 

Cabc — 7 abc ^ a b C (62) 

in case of a frame base. Here, 7 are the frame coefficients (including torsion), the symbol n denotes the contortion 
coefficients, while v denotes the nonmetricity coefficients. 

In the current version of ITENSOR, all these coefficients have now been implemented. While it is recognized 
not to be mathematically exhaustive, as a test the following relationships can be verified. First, when torsion is 
introduced: 


f-,ij - f;ji = - Kji k ) / )fc = 

if r is antisymmetric in its covariant indices. For nonmetricity, the following applies: 

Qij\k — 


(63) 

(64) 


where g denotes the metric tensor. 

3.1.6 Contraction of non-symmetric tensors 

Previously, when the ITENSOR package contracted tensor expressions, it did not take into account index ordering. 
Thus, for instance, e ab Tb c and e 6 a Xf, c both yielded the same result, even if e was declared as an antisymmetric tensor. 
The new implementation relies on the canform function to determine the canonical ordering of indices and any 
resulting sign changes before executing the contraction, thereby ensuring that the contraction operation provides 
consistent results. 

3.1.7 Exterior algebra 

Thanks in a very large part to the work of Valery Pipin, ITENSOR now has capabilities relating to fully antisym¬ 
metric covariant tensors. Paralleling the functionality of Maxima’s cartan package (which introduces an index-free 
formalism for differential forms), ITENSOR can now compute wedge products, contraction of an exterior form with 
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a vector, and the exterior derivative. Furthermore, ITENSOR also offers the capability to compute the Lie-derivative 
of a tensor of arbitrary valence. 

Since ITENSOR is a tensor algebra package, it appears natural to use the “tensorial” definition (TjTj!) of exterior 
forms. Many problems, however, can be expressed more naturally when the “geometric” definition S3 is employed. 
For this reason, a flag has been implemented in ITENSOR, which allows the user to use whichever definition is more 
applicable to the problem at hand. The igeowedgejflag variable is set to false by default, but if it is set to true, 
the “geometric” interpration is employed. 

3.1.8 Other changes 

One other major change that was implemented in the ITENSOR package was the consistent renaming of all objects 
and functions in order to conform to the naming conventions of commercial MACSYMA. The new names also make 
it a great deal less likely for ITENSOR names to conflict with names in the core Maxima system, or in other Maxima 
packages. 

3.2 Component tensor manipulation 

The component tensor manipulation package of Maxima, CTENSOR, views tensors as matrices. The main purpose 
of this package is to compute frequently used tensors in the coordinate system of the user’s choice, given the metric 
that the user provides. As an example, CTENSOR can be used to compute the Ricci-tensor of a specific metric and 
show, after simplification, whether or not the metric is that of an empty spacetime. 

Until recently, the CTENSOR package was largely non-functional. If the user tried to follow simple examples and 
invoked the csetup function, the program became apparently non-responsive, and Maxima had to be shut down. 
This problem turned out to be entirely superficial in nature: other components of Maxima failed to flush output 
buffers, causing an apparent lockup of the system when Maxima was waiting for user input, but the corresponding 
prompts were not yet displayed. 

Most of the work in CTENSOR since involved new features, most notably the support for frames, torsion, 
nonmetricity, and the algebraic classification of metrics. 

The CTENSOR package now has the capability to compute all the important tensors used in differential geometry. 
These include the Christoffel-symbols of the first and second kind, the Riemann tensor, the Ricci-tensor, the Einstein 
tensor, and the Weyl tensor. These tensors are calculated after the user has specified the coordinate variables, and 
entered the metric tensor in the form of a matrix, with elements that are functions of the coordinates. Special 
simplifications are applied when the metric tensor is diagonal, speeding up calculations considerably. These tools 
make many useful calculations possible; for instance, one demonstration file shows how Schwarzschild’s spherically 
symmetric solution to Einstein’s field equations can be derived. 

3.2.1 Using a frame base 

A rigid frame is characterized by its basis vectors e( Q ), where labels in parentheses are basis vector labels, not 
coordinate indices. We can, of course, also express the coordinates of the basis vectors relative to some coordinate 
system: the components of these (contravariant) vectors will be e^ a y 

In its current version, the behavior of the CTENSOR package is now controlled by a flag, cframe_f lag, which 
determines whether the Christoffel-symbols or the connection coefficients are used when computing curvature. 

When frames are used, the user specifies the metric by providing the covariant components of the frame base. 
(This is consistent with the convention of specifying a metric in a coordinate base by the covariant components of the 
metric tensor.) The user must also enter the frame metric. A function is provided that computes the contravariant 
frame base components, and also calculates the metric tensor. 

Note that no inverse functionality is provided: CTENSOR cannot compute a frame base if the metric is given. 
Though in principle it is possible to establish a set of equations that an orthonormal frame base must satisfy, and 
from these compute a frame base, this blind, “brute force” method rarely, if ever, results in a frame base that is of 
any use. 
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3.2.2 Algebraic classification of metrics 

New functions have been added to CTENSOR that can compute the covariant and contravariant components of the 
Newman-Penrose tetrad and, using this tetrad and the Weyl-tensor, compute the Petrov-classification of metrics. 
The Petrov algorithm is a direct implementation of the algorithm described in JHJ . This package has so far been 
tested with some simple cases, and it correctly obtained the Petrov type of the Schwarzschild, Kerr and anti-deSitter 
metrics, Godel’s Universe, the Taub-Nut plane vacuum, and Allnutt’s perfect fluid metric. Some of these results 
were cross-checked against results obtained from the dedicated computer algebra system SHEEP, used for research in 
general relativity. That said, it should be emphasized that to date, the petrov () function has only been tested in 
these simplest of cases and a comprehensive test suite that can exercise all branches in the decision tree used by this 
function remains to be written. 

The actual algorithm for the Petrov classification is implemented by the following subroutine, written in the 
Maxima language: 

petrovO :=block 

( 

[ 


0 , 

’N, 

’II,’ 

III, 

’D, 

MI, 

1—1 
1—1 

^1 

MI, 

’I, 

’ I , 

11, 

MI, 

13, 

14,15, 

’N, 

’I, 

’ I , 

19, 

MI, 

21, 

13,23, 

III,: 

19, 

11, 

27, 

7, 

23, 

15,31] 


P: 1 

], 

if psi[4]#0 then P:P+1, 
if psi[3]#0 then P:P+2, 
if psi[2]#0 then P:P+4, 
if psi[l]#0 then P:P+8, 
if psi[0]#0 then P:P+16, 

if numberp(T [P]) and T[P]>0 then 

( 

if T [P]=7 then 

if ratsimp(psi[3]~2-3*psi[2]*psi[4])=0 then ’D 
else ’II 

else if T[P]=11 then 

if ratsimp(27*psi [4]~2*psi [1]+64*psi[3]~3)=0 then ’II 
else ’I 

else if T[P]=13 then 

if ratsimp(psi[1]~2*psi [4]+2*psi [2]~3)=0 then ’II 
else ’I 

else if T[P]=14 then 

if ratsimp(9*psi[2]~2-16*psi[1]*psi[3])=0 then ’II 
else ’I 

else if T[P]=15 then 

if ratsimp(3*psi [2]~2-4*psi [1]*psi[3])=0 and 

ratsimp(psi[2]*psi [3]-3*psi[1]*psi[4])=0 then ’II 
else ’I 

else if T[P]=19 then 

if ratsimp(psi[0]*psi[4]~3-27*psi [3]~4)=0 then ’II 
else ’I 

else if T[P]=21 then 

if ratsimp(9*psi [2]~2-psi[4]~2) then ’D 
else ’I 

else if T[P]=23 then block 
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( 


[I:ratsimp(psi[0]*psi[4]+3*psi [2]~2)], 

if 1=0 and ratsimp(4*psi[2]*psi [4]-3*psi[3]~2)=0 then ’III 
else block 
( 

[J:ratsimp(4*psi[2]*psi[4]-3*psi [3]~2)], 

if ratsimp(psi[4]*I~2-3*J*(psi[0]*J-2*psi[2]*1))=0 then ’II 
else ’I 

) 

) 

else if T[P]=27 then 

if ratsimp(psi[0]*psi[3]~2-psi[1]~2*psi[4])=0 then 
if ratsimp(psi[0]*psi [4]+2*psi[1]*psi[3])=0 then ’D 
else if ratsimp(psi[0]*psi[4]-16*psi[1]*psi[3])=0 then ’II 
else ’I 
else block 
( 

[I:ratsimp(psi[0]*psi[4]+2*psi [1]*psi [3])], 
if 1=0 then block 
( 

[J:ratsimp(-psi [0]*psi[3]~2-psi[1]~2*psi[4])], 
if J=0 then ’III 

else if ratsimp(I~3-27*J~2)=0 then ’II 
else ’I 

) 

else ’I 

) 

else block 

( 

[H:ratsimp(psi[0]*psi[2]-psi[1] ~2)] , 
if H=0 then 

if ratsimp(psi[0]*psi [3]-psi [1]*psi [2])=0 then 
if ratsimp(psi[0] *psi [4]-psi[2]~2)=0 then ’N 
else ’I 
else block 
( 

[E:ratsimp(psi[0]*psi[4]-psi[2] ~2) ] , 
if E=0 then 

if ratsimp(37*psi[2]~2+27*psi[1]*psi [3])=0 then ’II 
else ’I 
else block 
( 

[ 

A:ratsimp(psi[1]*psi [3]+psi [2] ~2) , 

I: ratsimp(E-4+A) 

], 

if I#0 and ratsimp(I~3-27*(psi[4]*H-psi[3]~2*psi[0]+ 

psi[1]*psi[2]*psi[3]+psi[2]*A)~2)=0 then ’II 

else ’I 

) 

) 

else block 

( 

[I:ratsimp(psi [0]*psi[4]-psi[2]~2-4*(psi[1]*psi[3]+psi[2] ~2 ))], 
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if 1=0 then 

if ratsimp(psi[4]*H-psi[3]~2*psi[0]+psi[1]*psi[2]*psi[3]+ 
psi[2]*(psi[1]*psi[3]+psi[2]~2))=0 then ’III 

else ’I 
else 

if ratsimp(psi[0]~2*psi [3]-psi [0]*psi [1]*psi[2]-2*psi[1]*H)=0 then 
if ratsimp(psi [0]~2*I-12*H~2)=0 then ’D 
else if ratsimp(psi[0]~2*I-3*H~2)=0 then ’II 
else ’I 
else block 
( 

[J:ratsimp(psi [4]*H-psi[3]~2*psi[0]+psi[1]*psi[2]*psi[3] + 
psi [2] * (psi [1] *psi [3] +psi [2] ~2) )] , 
if J#0 and ratsimp(I~3-27*J~2)=0 then ’II 
else ’I 

) 

) 

) 

) 

else T[P] 

); 


For the simpler cases, the metric type is obtained directly from the lookup table T. For the more complicated 
cases, additional coefficients are computed on an “as and when needed” basis, in order to minimize the amount of 
computation required. 

3.2.3 Torsion and nonmetricity 

Torsion and nonmetricity can be represented as additional components in the connection coefficients. As defined 
in <ED and (PI). the connection coefficients can incorporate torsion and nonmetricity using the contortion and 
nonmetricity coefficients. Two functions have been added to the CTENSOR package to compute these coefficients 
from a user-supplied torsion tensor and nonmetricity vector. 

3.2.4 Predefined metrics 

Another useful feature found in the commercial MACSYMA version of CTENSOR was predefined metrics. Through 
the ct_coordsys () function, CTENSOR had the ability to set up a variety of commonly used metrics in a coordinate 
or frame base. 

This capability has also been reproduced in the current version of CTENSOR. Specifically, the following metrics 
have been implemented: 

Name/Coordinates Metric (gr^/Frame (e( a p) 

cartesian2d (i o\ f 1 o\ 

x,y ^0 lj’ ^0 1 J 


polar 

r,<f> 


f 1 0 \ / cos <j> 

^ 0 r / ’ \ s i n 4 1 


r cos (f> 


f e 2 (cosh 2 u — cos 2 v) 

\ 0 


0 

e 2 (cosh 2 u — cos 2 v) 


( e sinh u cos v 
e cosh u sin v 


—e cosh u sin v 
e sinh u cos v 


elliptic 

u, v 

e is constant 
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confocalelliptic 
u, v 

e is constant 



0 



ev 

eu(l — v 2 ) 


eu 

ev(l — u 2 ) 


V (u 2 -l)(l-v 2 ) \/ ( u 2 -l)(l-v 2 ) 



( 

° 2 ) 


I (cosh —cos u) z 

bipolar 

\ 0 

e / 

(cosh d — cos u) 2 / 

u, V 



e is constant 

/ — e sin it sinh v 

e(l — cos u cosh v) \ 


( (cosh i;—cos u) 2 

(cosh v—cos u) 2 \ 


\ e(cos u cosh v — 1) 

— e sin u sinh v / 


\ (cosh v — cos u) 2 

(coshv-cos / 


parabolic 
u, v 


2 i 2 

U + V 


0 

it 2 + V 2 


— 


cartesian3d 

x,y,z 


/I 0 0\ /I 0 0\ 

0 1 0,0 1 0 

\0 0 1 / \0 0 1 / 


polarcylindrical 
r, 9 ,2 


/I 0 0 \ / cos 6* — r sin 0 0 \ 
10 r 2 0 I , j sin 8 r cos 6 0 ] 
\0 0 1 / \ 0 0 1 / 


ellipticcylindrical 
u, v, z 

e is constant 


/ e 2 (sin 2 v + sinh 2 u) 0 0 \ 

j 0 e 2 (sin 2 v + sinh 2 u) 0 I 

Vo 0 1/ 


( e sinh u cos v —e cosh u sin v 0 ' 
e cosh u sin v e sinh u cos v 0 

0 0 1 , 


confocalellipsoidal 
u, v, w 

e, f, g are constants 


(v — u)(w — u) 

4 (e 2 -u)(u-f 2 )(v,-g 2 ) 

0 

0 


(v — u)(w — v) 

4 (v-e' 2 )(v-f 2 )(v-g' 2 ) 


0 

0 

(w — u)(w — v) 

4(e‘ J —w)(w — f 2 )(w — g 2 ) 


(- 


(v — e 2 )(w — e 2 ) 

4( f 2 -e 2 )(g 2 -e 2 )(e 2 -u) 

(v-f 2 )(w-f 2 ) 

4 (f 2 -e 2 )(9 2 -f 2 )(u-f 2 ) 


w 


( V -g2 Kw - g 2 ) 

Ma' 2 -e 2 )(g 2 -f 2 )(g 2 -u) 


(u — e 2 )(w — e 2 ) 
4(f 2 —e 2 )(g 2 —e 2 )(e 2 —v) 

(u-f 2 )(w-f 2 ) 
Hf 2 -e 2 ){g 2 -f 2 ){v-f 2 ) 
( u g 2 )( w g2) 

4 (g' 2 -e 2 )(g 2 -f 2 )(g 2 -v) 


(u — e 2 )( v — e 2 ) 
4{f 2 -e 2 )(g 2 -e 2 ){e 2 -w) 

(u-/ 2 )Q-/ 2 i 

Mf 2 -e 2 )(g 2 -f 2 )(w-f 2 ) 


\ 


t (u-g 2 )(v-g 2 ) j 

"V 4 (g 2 -e 2 )(g 2 -f 2 )(g 2 -w) ) 
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bipolarcylindrical 

u , v, z 

e is constant 


(cosh v —cos u) 2 
2e 2 sin u sinh i>(l — cos u cosh v) 
(cosh v —cos u )-“* 

0 


2e 2 sin u sinh v(l — cos u cosh v ) 
(cosh v — cos u)2 
e 2 

(cosh v — cos u 

0 


0 

0 

e 2 

(cosh v — cos u) 4 


e sin u sinh v 
(cosh v —cos 
e (cos u cosh i; — 1) 
(cosh v —cos u) 2 

0 


e(l—cos u cosh v ) 
(cosh t) — cos 
— e sin u sinh v 
(cosh v — cos u ) a 

0 


0 

0 

_e_ 

(cosh t)—cos 


paraboliccylindrical 
ti, v, z 


/ u 2 + v 2 0 0\ / u —v 0\ 

I 0 u 2 + v 2 0 ), I v u 0 I 

\ 0 0 1 / \0 0 1 / 


paraboloidal 

u, v, (p 


/mV 0 0 \ 

0 m 2 + v 2 0 

\ 0 0 m 2 +v 2 J 


—UV Sill (p 
UV COS <p 

0 


V COS (p 
v sin <p 
u 



conical 

M, V, W 

e, f are constants 


( (v 2 — u 2 )w 2 
(u 2 — e 2 ) (u 2 — f 2 ) 

0 

0 


0 

0 

1 


/ 

V 


vw 

e f _ 

— uwy/ 

ei/(/ 2 - d)(» 2 - e 2 ) 

UW^J V 2 — f 2 

e^(f 2 -e 2 )(u 2 -f 2 ) 


ef _ 

— vwy/e^—v? 
e^/(f 2 - e 2 )(v 2 - e 2 ) 
vwyj u 2 — f 2 

f\J (f 2 -e 2 )(v 2 -f 2 ) 


' f _ 

(u 2 — e 2 )(v 2 — e 2 ) 
e\/e 2 -f 2 
y/(u 2 -f 2 )(v 2 -f 2 ) 


\ 


toroidal 

u, v, cp 

e is constant 


C sinli V 
(cosh v — cos u)2 

0 

0 


— e sin 4> sinh v 
(cosh v—cos u ) 
e cos </) sinh v 
(cosh v—cos u) 

0 


^cosh”iJ—cos”uJ^ 

0 


— e cos <f) sin u sinh v 
(cosh v — cos u)'* 

— e sin <fi sin u sinh('u) 

(cosh v — cos it) 2 
e(cos u cosh v — 1) 
(cosh v — cos u 


0 

0 

e 2 

(cosh v —cos u) 2 . 

-e cos <j)(cos u cosh v — 1) 
(cosh v —cos u)^ 

-e sin (f>(cos u cosh v — 1 ) 
(cosh v —cos u)^ 

— e sin u sinh v 
~(cosh~v--cosu)^ 


spherical 

r, 6, p 


/10 0 \ /I 0 0 \ 

Or 2 0 , 0 r 0 

\0 0 r 2 sin 2 6 J \0 0 rsinfl/ 


oblatespheroidal 

u, v, (p 

e is constant 


( e 2 (sin 2 v + sinli 2 m) 

s_ 

( |e| \J sin 2 v + sinh 2 u 


0 0 

e 2 (sin 2 v + sinh 2 u) 0 

0 e 2 cosh 2 m cos 2 v, 

0 _ 0 

|e| \/ sin 2 v + sinh 2 u 0 

0 lei cosh m| cos 
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oblatesplieroidalsqrt 

u, v, <f> 

e is constant 


prolatespheroidal 

u, v, <j> 

e is constant 


‘(u 2 -v 


0 

2 (u 2 -v 2 ) 

u' 2 -l 

0 


^ |e| \Ju 2 —v 

\J u 2 — 1 
0 
0 


0 0 \ 

e| \J u 2 —v 2 n 

0 \euv\) 


/ e 2 (sin 2 v + sinh 2 it) 0 

( 0 e 2 (sin 2 v + sinh 2 u ) 

\0 0 


S ) 

e 2 sin 2 v sinh 2 u ) 


( |e| \J sin 2 v + sinli 2 u 

l 


0 _ 

|e| \J sin 2 v + sinli 2 it 
0 


: ) 

|esinhitsint;| / 


prolatespheroidalsqrt 

u, v, <j> 

e is constant 


e 2 (v 2 —u 2 ) 
1 -u 2 


0 

2 (v 2 — u 


0 


0 
0 

e 2 (l - u 2 )(v 2 - 1), 


/ [e| \/v 2, —u 2 

o 


V 0 


0 


|e| yj v 2 —u 2 
\/v 2 — 1 
0 


o > 

0 

\e\y/(l-u>)(v*-l)J 


/ (a 2 cos 2 (f) + b 2 sin 2 (/>) sin 2 6+ 

( a 2 cos 2 <j> + b 2 sin 2 <j> — c 2 )- 

(b 2 

— a 2 ) cos 0sim(> • \ 

ellipsoidal 

r, e, <j) 

+c 2 cos 2 6 

■r cos 9 sin 9 


■r sin 2 9 

( a 2 cos 2 <j> + b 2 sin 2 <j> — c 2 )- 

r 2 ((a 2 cos 2 (j>+ b 2 sin 2 <j>) cos 2 9+ 

(b 2 

— a 2 ) cos 4> sin <f> ■ 

■r cos 9 sin 9 

+c 2 sin 2 9) 


■r 2 cos 9 sin 9 

a, b , c are constants 

( b 2 — a 2 ) cos (j> sin (j>- 
\ -r sin 2 9 

( b 2 — a 2 ) cos (f> sin </>■ 

■r 2 cos 9 sin 9 

(a 2 sin 2 <j> + b 2 cos 2 <j>) ■ 1 
■r 2 sin 2 9 J 


no frames 



/I 

0 

0 



/I 

0 

0 

°\ 

cartesian4d 

0 

1 

0 

0 


0 

1 

0 

0 

x,y,z,t 

0 

0 

1 

0 


0 

0 

1 

0 


Vo 

0 

0 

l) 


Vo 

0 

0 

1/ 



/I 

0 

0 

° ^ 


/I 

0 

0 

° \ 

spherical4d 

0 

r 2 

0 

0 


0 

r 

0 

0 

r,9,r),<t> 

0 

Vo 

0 

0 

r 2 sin 2 6i 

0 

0 

2-2 • 2 n i 

r sin 77 sin 6 / 

5 

0 

Vo 

0 

0 

r sin 6 1 

0 

0 

r sin rj sin 9 ] 


exteriorschwarzschild 

/ 

2m—r 

0 

0 

0 V 



0 

0 

0 

t, r , 9, (f> 


0 

r 

r—2m 

0 

0 


0 

\/ r-2m 

0 

0 

m is constant 


0 

0 

r 2 

0 

' 

0 

0 

r 

0 

r > 2m 

V 

0 

0 

0 

r 2 sin 2 9 / 


V 0 

0 

0 

r sin 9 


interiorschwarzschild 
t, z, u, v 
m is constant 
t < 2m 


( 2m —t 

0 

0 

0 

2m — t 

0 

0 

t 2 

0 ^ 
0 

0 


^ Y 2m—t 

0 

0 

0 \ 

t 

0 

5 

0 

0 

/ 2m— t 

V t 

0 

0 

f 

0 

0 

V 0 

0 

0 

t 2 sin 2 u J 


V 0 

0 

0 

t sin u J 
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2 mr — r 2 — q. cos 2 0 
r 2 + a z cos^ 0 

0 


0 


— 2amr sin 2 6 
r z -\-a z cos'* # 

0 


kerr_newman 

t,r,9,<t> 

a, m are constants 



The signature of all these metrics is (+,+) with the exception of the Schwarzschild and Kerr-Newman metrics 
that have a (—, +,...+) signature. The ct_coordsys function also has the ability to add additional “flat” dimensions 
with either a Minkowski or a Euclidean signature. 

3.3 Algebraic tensor manipulation 

In addition to CTENSOR and ITENSOR, the commercial MACSYMA system included a third tensor package, 
ATENSOR. The purpose of this package was to capture the algebraic properties of “general tensor algebras”. This 
package allowed the user to define a tensor algebra by specifying the algebra’s commutation properties. A variety of 
algebras could be explored this way, including Clifford algebras, Grassmann algebras, and Lie algebras. 

To make a long story short, unlike CTENSOR and ITENSOR which, while in various states of disrepair, were 
nevertheless present in GPL Maxima, the ATENSOR package was completely missing. This has now been corrected, 
and an initial implementation of ATENSOR, fully conformant to the documentation of commercial MACSYMA, is 
now included with GPL Maxima. 

3.3.1 Commutators and abstract algebra types 

The ATENSOR package defines different algebra types based on their (anti)commutators. Several algebra types are 
recognized: 

Algebra type Commutator 

Universal No commutation rules 

Grassmann u ■ v + v ■ u = 0 

Clifford u ■ v + v ■ u = 2 f s {u, v) 

Symmetric u • v — v ■ u = 0 

Symplectic u ■ v — v ■ u = 2f a (u, v) 

Lie enveloping u ■ v — v ■ u = 2v a (it, v) 

Here, f s stands for a scalar valued symmetric function (i.e., f s (u,v) = f s (v,u)), f a stands for an antisymmetric 
scalar valued function ( f a {u,v ) = —f a {v,u)) and v a is an antisymmetric vector valued function. 

The actual values of the functions / s , f a , and v a are determined by the algebra’s type and dimensionality. The 
ATENSOR package uses the matrix af orm in which the function values for base vector arguments are stored. For the 
vector-valued function v a , the matrix values are taken to be base vector indices. When the algebra type is selected 
by the user, this matrix is preinitialized appropriately. 

For this purpose, when a Clifford, Symplectic, or Lie enveloping algebra is selected, the user can optionally enter 
values that determine the algebra’s dimensionality. For a Clifford algebra, up to three values can be used, specifying 
the positive, degenerate, and negative dimensions of the algebra, respectively. For a symplectic algebra, the number 
of regular and degenerate dimensions can be specified; for Lie enveloping algebras, a single number, the algebra’s 
dimensionality, can be entered. 
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For a Clifford algebra, the matrix af orm is initialized as diag(l... 1, 0 ... 0, —1... — 1), where the number of l’s, 
0’s, and —l’s correspond with the number of positive, degenerate, and negative dimensions, respectively. 

For a symplectic algebra of (regular) dimension n, aform is initialized to an n x n matrix with its off-diagonal 
values set to 1 or —1 depending on whether they represent an odd or even index permutation. If the algebra also 
has degenerate dimenions, the appropriate number of null columns and rows are added. 

Lastly, for Lie enveloping algebras aform is initialized to an antisymmetric matrix whose elements are defined as 
ciij = (mod„(2n + 2 — * — j) + 1)perm(i, j) where perm(«, j) is a permutator function that gives +1 if — 

1, i + 1,... j — 1, j + 1,... n) is an even permutation of the sequence (1... n), and —1 otherwise. For instance, the 
initialization function call init_atensor(lie_envelop,3) produces the following matrix: 

0 3-2 

-3 0 1 

2-10 


It is known (see, for instance, 
with the algebra of quaternions, 
table as 


P) that the Clifford algebra of 0 positive and —2 negative dimensions corresponds 
The ATENSOR package can correctly reproduce the quaternionic multiplication 


1 

Vl 

V2 

Vl ■ Vl 

Vi 

-1 

Vl ■ Vi 

- Vl 

Vi 

-Vl ■ v 2 

-1 

Vl 

VI ■ V2 

Vi 

—Vl 

-1 


where the base vectors v±, tq, and their product iq • tq serve as the three quaternionic imaginary units. This was 
one of the tests used to ascertain that the ATENSOR package produces mathematically correct results. 


4 Conclusions 

The Maxima tensor packages remain a work in progress. It is the author’s desire to maintain these packages in good 
working order so long as GPL Maxima itself remains actively maintained. It is hoped that the tensor packages will 
again be accepted by a user community, and as they are actively “field tested”, it will be possible to make them 
more robust and mathematically accurate. 
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